% Latest series of computation for ocean stuff
% Aim is to read and understand Alexei's file format
% of source file
%further aims to compute voltages - project the electric field vectors etc
%akola 19.11.2005

%fid = fopen('e:\manoj\ocean\x3dg\1x1_winter_2002.source','rt');
%fid = fopen('c:\manoj\ocean\x3dg\1x1_summer_2002.source','rt');
%fid = fopen('c:\manoj\ocean\source\180x360\monthly\apr00.in.source','rt');
%fid = fopen('e:\manoj\ocean\x3dg\temp.source','rt');
fid = fopen('E:\manoj\ocean\Source\180x360\monthly\dec98.source','rt');
%Caution Alexei used the wrong transport values (1/3rd) for this
%Source file...... 27.9.5
flag = 1;

while flag == 1,
    A = fgetl(fid);
    if length(findstr(A,'ARRAY')) == 1 & findstr(A,'ARRAY') == 5,
        flag = 0;
    end;
end;

[D1,count] = fscanf(fid,'%f',[360,180]);

flag = 1;

while flag == 1,
    A1 = fgetl(fid);
    if length(findstr(A1,'ARRAY')) == 1 & findstr(A1,'ARRAY') == 5,
        flag = 0;
    end;
end;

[D2,count] = fscanf(fid,'%f',[360,180]);

flag = 1;
while flag == 1,
    A1 = fgetl(fid);
    if length(findstr(A1,'ARRAY')) == 1 & findstr(A1,'ARRAY') == 5,
        flag = 0;
    end;
end;

[D2,count] = fscanf(fid,'%f',[360,180]);

fclose(fid);

% old U V fields from the file sent Kuvsh
% data = load('c:\manoj\ocean\velocities\uV_World_2002.txt');
% data = load('c:\manoj\ocean\velocities\UV_World_summer_new_2002.txt');
data = load('c:\manoj\ocean\velocities\UV_World_summer_wrong_2002.txt');
U1 = reshape(data(:,3),[360,160]);
V1 = reshape(data(:,4),[360,160]);
L1 = V1 == 9999;
L2 = U1 == 9999;
L = L1|L2;
U1(L2) = 0;
V1(L1) = 0;
UU = [zeros([360,10]) U1 zeros([360,10])];
VV = [zeros([360,10]) V1 zeros([360,10])];

UU = UU';
VV = VV';
load c:\manoj\temp\bz.mat;
Bz = flipud(Bz');

load c:\manoj\ocean\velocities\world_2002.jex;
U = world_2002(:,3);
V = world_2002(:,4);
U = reshape(U,[360,180]);
V = reshape(V,[360,180]);

%---------------------------------------------------------

% Define Grids
% The original grids were assumed to be 1x1 degree (center)
% change them if NOT
    
worldmap('lo',[0 50],[110 180]);
h=textm(13.6,144.9,' Guam');
set(h,'Color','b','FontSize',14);
h=textm(15.8,121.6,' Philippines');
set(h,'Color','b','FontSize',14);

plotm(15.8,121.6,'r.','MarkerSize',20);
plotm(13.6,144.9,'r.','MarkerSize',20);
h=textm(28.2,177.4,' Midway');
set(h,'Color','b','FontSize',14);

plotm(28.2,177.4,'r.','MarkerSize',20);
h=textm(26.1,127.8,' Okinawa');
set(h,'Color','b','FontSize',14);

plotm(26.1,127.8,'r.','MarkerSize',20);
h=textm(32.5,130.1,' Amakusa');
set(h,'Color','b','FontSize',14);

plotm(32.5,130.1,'r.','MarkerSize',20);
h=textm(30.2,128.2,' Shanghai'); 
set(h,'Color','b','FontSize',14);

plotm(35.3,139.3,'r.','MarkerSize',20);
h=textm(35.3,139.3,' Ninomiya'); 
set(h,'Color','b','FontSize',14);

plotm(30.2,128.2,'r.','MarkerSize',20);

linem([15.8,13.6],[121.6,144.9],'r.-');
%--------------------------------------------------------
% Project the vector data to a line find scalar lengths & integrate - to
% find the voltage 19.11.2005 akola

load c:\manoj\projects\cable\5aug06\ecco1x1_comp_fields.mat Et Ep; % The simulated Electric field from ECCO 1992-2002 fields
load c:\manoj\projects\cable\5aug06\IUIV199.mat timeaxU;
TimeID = TimeaxU;

%Guam-Philipines
% clat = [15.8,13.6]; % 
% clon = [121.6,144.9];

%Guam - Okinawa
% clat = [26.1,13.6]; % 
% clon = [127.8,144.9];

%Guam - Midway
% clat = [13.6,28.2];
% clon = [144.9,177.4];
% clat = [28.2,13.6];
% clon = [177.4,144.9];

%Ninomiya-Okinawa
clat = [35.3,26.1];
clon=[139.3,127.8];


dateE = TimeID/(3600*24) + datenum(1992,1,1,6,0,0); % 1992 Jan 1, 6 AM is the starting point

for i = 1:132,
    [Vp,r1,lat,lon] = mapprofile(flipud(squeeze(Epr(i,:,:)))*1e6,[1,90,0],clat,clon);
    [Vt,r1,lat,lon] = mapprofile(flipud(squeeze(Etr(i,:,:)))*1e6,[1,90,0],clat,clon);
    pvec = [-diff(clat),-diff(clon)];
    pvec = pvec/(sqrt(sum(pvec.^2)));
    pvec=repmat(pvec',[1,length(Vp)])';
    %dvec = [-Vt Vp];
    dvec = [Vt Vp];
    sc =dot(pvec',dvec');
Vol(i) = trapz(deg2km(r1),sc);
end;



%[map0,maplegend0] = maptrims(flipud(circshift(squeeze(IV(1,:,:)),[0,180])),[1,90,0],[0 50],[110 180]);
%worldmap(map0,maplegend0);

[map0,maplegend0] = maptrims(flipud(squeeze(Etr(i,:,:)))*1e6,[1,90,0],[0 50],[110 180]);
worldmap(map0,maplegend0);
